#include "ms_sphere.h"

using namespace std;

/* ******************************************************************* */
sphere::sphere( const double r )
{
	X[0] = X[1] = X[2]	= 0.;
	R		= r;
}
/* ******************************************************************* */
sphere::sphere( const double x[], const double r )
{
	X[0]	= x[0];
	X[1]	= x[1];
	X[2]	= x[2];
	R		= r;
}
/* ******************************************************************* */
double	sphere::Distance(const sphere & s)  const {
	const double d =	sqrt(	(s.X[0]-X[0])*(s.X[0]-X[0])
							+	(s.X[1]-X[1])*(s.X[1]-X[1])
							+	(s.X[2]-X[2])*(s.X[2]-X[2])	);
	if( d - s.R - R <= 0. ) return -1.;
	else return 0.;
}
/* ******************************************************************* */
bool	sphere::Collision(const sphere & s) const {
	return (Distance(s) < 0.);
}
/* ******************************************************************* */
double	sphere::IntersectionVolume(const sphere & s) const {
	const double d =	sqrt(	(s.X[0]-X[0])*(s.X[0]-X[0])
							+	(s.X[1]-X[1])*(s.X[1]-X[1])
							+	(s.X[2]-X[2])*(s.X[2]-X[2])	);
	const double	RR = ( ( R < s.R ) ? s.R : R ),		/* larger of the two radii */
					rr = s.R + R - RR;					/* smaller of the two radii */
	if( d>s.R+R ) return 0.;
	const double		PI1		= acos(-1.),
						PI3		= PI1/3.,
						PI43	= 4.*PI1/3.;

	if ( d + rr <= 1.0001*RR ) return PI43*rr*rr*rr; /* small sphere in big sphere */
	if ( d*d+rr*rr <= RR*RR ) {
		const double	a		= 0.5 * ( (RR*RR-rr*rr+d*d)/d );
		const double	b		= d-a;
		return PI3*( (RR-a)*(RR-a)*(2.0*RR+a) + (rr-b)*(rr-b)*(2.0*rr+b) );
	}
	const double	a		= 0.5 * ( (RR*RR-rr*rr-d*d)/d );
	const double	b		= d+a;
	return PI3*( (rr*rr-a*a)*(2.0*rr-a) + (RR-b)*(RR-b)*(2.0*RR+b) );
}
/* ******************************************************************* */
void sphere::DrawVoxelImage(	VoxelImage & vx_img,
								const unsigned char col ) const 
{
	int iz0 = vx_img.GetVoxelCoZ( GetZMin() ),
		iz1 = vx_img.GetVoxelCoZ( GetZMax() );
	if( iz1 < 0 ) return;
	if( iz0 >= vx_img.Nz) return;
	if( iz0 < 0 ) iz0 = 0;
	if( iz1 >= vx_img.Nz) iz1 = vx_img.Nz-1;
	double z = vx_img.GetZ(iz0); // (double)iz0*vx_img.h + 0.5*vx_img.h - 0.5*vx_img.wx;
	double r_xy ;
	for(int iz=iz0; iz <= iz1; iz++)
	{
		r_xy = R*R - (z-X[2])*(z-X[2]);
		if(r_xy > 0.) {
			r_xy = sqrt(r_xy);
			if(vx_img.DrawCircleZ( (double[]){X[0],X[1],z}, r_xy, col )  != 0);
				// printf("ERROR in slice %d of %d %16.8lf !\n", iz, iz1, r_xy);
		}

		z += vx_img.h;
	}
};
